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We study head-on collisions of boson stars in three dimensions. We consider evolutions of two 
boson stars which may differ in their phase or have opposite frequencies but are otherwise identical. 
Our studies show that these phase differences result in different late time behavior and gravitational 
wave output. 
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I. INTRODUCTION 

Models of compact "stars" described by non-fluid 
sources serve both as probes of strong field gravity as 
well as exotic alternatives to the standard "regular" fluid 
stars. Among these non-standard star models, a geon -a 
self gravitating star consisting of electromagnetic fields- 
was first considered by Wheeler [lj . Despite its appeal, no 
stable geon could be constructed and so interest dimin- 
ished. However, borrowing from this idea, Kaup was able 
to construct the first stable configuration -referred to 
as a Klein-Gordon geon- by considering a U(l) classical 
scalar field @- Soon after Kaup's construction, Ruffini 
and Bonazzola revisited the solution, adopting a quan- 
tum point of view for the matter fields though maintain- 
ing the geometry as a classical entity @. By employing 
a Hartree-Fock (multi-particle) approximation, they re- 
obtain Kaup's equations (and hence solutions) providing 
an alternative view to the resulting stars. These solutions 
are currently known as boson stars and describe a family 
of self-gravitating scalar field configurations within gen- 
eral relativity. 

Boson stars have been studied in many different con- 
texts (see [UjH] for reviews). They have been proposed 
as models for dark matter (see [1] for a recent review) and 
more recently have been used as convenient models for 
compact objects including black holes Q. One can con- 
sider strongly gravitating compact-star spacetimes with- 
out the worries associated with regular matter -such as 
shock fronts and discontinuities in the fluid variables- 
making boson stars very useful probes of strong-field gen- 
eral relativity. 

Given this utility, a number of efforts have modeled 
the dynamical evolution of boson stars. For instance, 
they have been considered in one and two dimensions 
for studies of critical phenomena @, 0] . In three dimen- 
sions, studies of the dynamics of perturbed boson stars 
have bee n p resented in [Tol . [Tlj and preliminary colli- 
sions in [12j. Works examining boson star collisions in 
two dimensions were carried out in Q while in three di- 
mensions, in addition to our own, another effort is under 
way [HI]. Another work which studies the collisions of 



solitonic solutions of the Klein-Gordon equations in the 
absence of gravity with an attractive potential. 

In this paper we investigate several boson star configu- 
rations paying particular attention to head-on collisions 
in 3D. Our aim is to survey some of the phenomenol- 
ogy that can be found in these scenarios. Additionally, 
since one can employ this particular physical system to 
describe strongly gravitating compact stars, one can use 
them to model compact systems in general relativity in 
their early stages where tidal effects are not too severe. 
On a more speculative level, one can also examine the 
gravitational wave signature produced so that searches 
by gravitational wave detectors might place constraints 
on their existence. 

The paper is organized as follows. In Section [II] we de- 
tail the formalism we use for the Einstein-Klein-Gordon 
system. Section [IIII describes the numerical implementa- 
tion of the governing equations. In Section IV we present 
our results. We conclude in Section [V] adding some final 
comments. The computation of the initial data for a 
single boson star is described in the Appendix [Al while 
the Appendix [B] is devoted to some (qualitative) physi- 
cal considerations in order to explain the results of the 
simulations. 



II. THE EINSTEIN-KLEIN-GORDON SYSTEM 

The dynamics of a complex scalar field in a curved 
spacetime is described by the following Lagrangian den- 
sity (adopting geometrical units, i.e. G = c = 1) 
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g ab d a 4>d b <t> + V |0| 



(1) 



Q-balls also has relevance here [14|. These Q-balls are 



Here R is the Ricci scalar, g a b is the spacetime metric, <f> 
is the scalar field, <f> its complex conjugate, and V(|0| 2 ) a 
potential depending only on \<p\ 2 . Throughout this paper 
Roman letters at the beginning of the alphabet a, 6, c, .. 
denote spacetime indices ranging from to 3, while let- 
ters near the middle i,j,k,.. range from 1 to 3, denoting 
spatial indices. This Lagrangian gives rise to the equa- 
tions determining the evolution of the metric (Einstein 
equations) and those governing the scalar field behavior 
(Klein-Gordon Equations) . 
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A. Einstein Equations 

The variation of the action associated with the La- 
grangian ([T]) with respect to the metric g ab leads to the 
Einstein equations 



R 

Rab — -^g a b = 8nT a b, 



(2) 



where R a b is the Ricci tensor and T ab is the stress-energy 
tensor given by 



T a b = ^ [d a <j)db<p + d a (f> d b <j)] 

- \g a b[g cd d c j>d d( t> + V(\<i>\ 2 )] 



(3) 



The Einstein equations form a system of 10 non-linear 
partial differential equations for the spacetime metric 
components g a b- A convenient way to express the equa- 
tions can be obtained by the identity 

Rab = —-^g cd d c ddg a b + d( a T b ^ —T ca bY c 

+g cd g ef (d e g ca d f g db - T ace T bdf ) . (4) 
where we have introduced the quantities 

T a = g bc T abc (5) 
Tabc = - (d b g ac + d c g ab - d a gbc) ■ (6) 

To this point we have considered arbitrary coordinates 
x a that label spacetime points. In order to define an ini- 
tial value problem for the Einstein equations, a foliation 
defined by the hypersurfaces with x° = t = const is intro- 
duced with normalized normal vector n a — — V a f/| | V a t| |, 
where V a is the covariant derivative associated to g ab . 
The harmonic formulation of Einstein equations exploits 
the fact that the coordinates x a can be chosen satisfy- 
ing the harmonic V C V c x a — [16] or their generalized 
version [HI El] 



V c V c x a = H a (t,x l ) 



(7) 



where the functions H a amount to prescribing coordi- 
nate conditions. These expressions can be combined to 
re-express Einstein equations p|) in their generalized har- 
monic form, as [l?], Ell 



T 



g c d cd g a b + d a H b + d b H a = -16 it yT ab - — g a b 

+2 T cab H c + 2 g cd g e f [ d e g ac d f g bd - T ace T bdf ] (8) 

Notice that the partial derivatives of H a do not belong 
to the principal part of the system since they are pre- 
scribed spacetime functions. The principal part of ((Sj) 
is a well posed system of decoupled wave equations for 
the 10 metric components. It is important to stress that 
equations ([5]) with the generalized harmonic condition 



([7]) constitute an overdetermined system since T a can be 
obtained either from the choice of coordinates ([7]) or from 
the evolution of the metric itself ijHJ). In the free-evolution 
approach, equation ([5]) is adopted to update the metric, 
and equation {7J is regarded as a constraint on the full 
system. Possible constraint violations can be measured 
by introducing a new four-vector [47| 



2 Z a = -T a -H a (t,x i 



(9) 



Clearly, the physical solutions will be those satisfying the 
algebraic condition Z a = throughout the spacetime; we 
will refer to the Z a quantities as the Z-constraints from 
now on. Note that added flexibility, useful for numeri- 
cal purposes, is attained by adding these constraints to 
the equations. In particular, since the addition of these 
Z-constraints to the evolution system can change the sta- 
bility of the solutions against perturbations off the con- 
straint surface one can exploit this fact to one's advan- 
tage. Indeed, suitable terms can be added to the equa- 
tions in order to construct an attractor for the physical 
solutions in certain spacetimes, in such a way that small 
Z-constraint violations will be damped during the evolu- 
tion Q ■ With these damping terms, equation ([5]) can 
be written as 

g cd d cd9ab + d a H b + d b Ha= (10) 
+2 T cab H c + 2 g cd g e f [ d e g ac d f g bd - T ace T bdf } 

-2 do [n a Z b + n b Z a - g a bn c Z c ] - 16 7T ^T ab - g ab 

where we have introduced <jq, which is a free parameter 
that controls the damping of the Z-constraints (0 . 



B. Klein-Gordon Equations 

As mentioned, the variation of the Lagrangian |T]) with 
respect to the scalar field <j> leads to the Klein-Gordon 
(KG) equations 



n<j> = 



dV 



(11) 



where the box □ = g ah V a V b stands here for the wave 
operator on a curved background. For concreteness, from 
now on we will consider the free field case where the 
potential takes the form 



(12) 



with m a parameter that can be identified with the bare 
mass of the field theory, although it has units of inverse 
length. For our simulations we have fixed its value to 
unity. The potential (fT2"|) leads to the so-called minibo- 
son stars, because achievable stable configurations have 
small masses. More general terms can be included, such 
as the A |(/>| 4 self-interaction term introduced in [2l|, lead- 
ing to heavier boson stars which have masses and sizes 
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more relevant to astrophysical applications when consid- 
ering their compactification ratio. In the present work 
we restrict to the A = case, deferring to a future work 
the study of the effect of this self-interacting term. 



Additionally, as with the Z-constraints, the mathemat- 
ical properties of the evolution system can be changed 
-when the constraints are not exactly satisfied- by judi- 
ciously adding the C-constraints to the right-hand-side 
of the equations. For instance, in [19( these constraints 
are incorporated in the equations in such a way that 



C. Reduction of the complete system to first order the physical solutions (i.e., those satisfying Cj = Cj 



The second order Einstein-Klein-Gordon (EKG) equa- 
tions (fT0|| - (|ll|) . can be employed to obtain the future evo- 
lution of the spacetime if {g a b- l 4>-i9ab,t,4>,t} are provided 
at an initial hypersurface. For a numerical implemen- 
tation, it is also useful to further reduce the system to 
fully first order to take advantage of several numerical 
techniques devised to exploit stability theorems for these 
kind of systems [13, [H, 0, [2f| . Consequently, these tech- 
niques provide a clean path to discretize the system in a 
robust manner 0, [U • Thus we begin by reducing 
these equations to first order form. 

The reduction in time is achieved by introducing new 
independent variables related to the time derivatives of 
the fields. Following a similar route as in [l9(, we define 



Qab 



dc9ab 



n 



dc 



(13) 



The evolution equations for Q a b and II are now given by 
the EKG equations (|10lll[) . while the evolution equations 
for g a b and <f> are simply given by their definition (|13|) . 

The reduction to first order in space is made by in- 
troducing new independent variables encoding the first 
space derivatives as 



Diab = dig a b , $i = di< 



(14) 



The equations of motion for these quantities are obtained 
by applying the time derivative to their definition ([T4]) 
and imposing that time and spatial derivatives commute. 
Notice that in this step, just as with the Z-constraints, 
one encounters an overdetermined system since spatial 
derivatives can be obtained either by deriving the fields 
or from the evolution equation of the new quantities (fl4|) . 
This fact can be re-interpreted as adding new constraints 
to the system 

Ciab = dig a b - D iab = , Ci = di<f) - $i = , (15) 



Ciab = Cij a b = 0) are an attractor in certain spacetimes. 
This means that small C-constraint violations will also 
be damped during the evolution. Wc follow the same 
approach here for the full first order reduction of the 
Einstein Klein-Gordon equations by writing them as: 

dtgab = P k D kab -a Q ab (17) 
d t Qab = P k d k Qab - n -•"<•;, />.„.,, - CTll3 k d k g a b (18) 

- a d a H b - a d b H a + 2 a T cab H c 

+ 2ag cd {f j D ica D jdb - Q ca Qdb - g ef T aC eT bdf ) 

- -n c n d Q cd Q ab - a -y v D iab Q jc n c + cri/3 fc D fcab 

- 16na (T ab - 9 -fr) 

- 2 a cr [n a Z b + n b Z a - g ab n c Z c ] 

d t Diab = (3 k d k D ia b - a d l Q ab + a a 1 d l g ab (19) 

+ ^n c n d D icd Qab + a j 3k n c D lJC D ka b - a aiD iab 

d t <j) = p k <P k -all (20) 
dtli = p k d k U-aY'd l ^ J ;-aiP k d k <t> (21) 

- n c n d Q cd - a f 3 $ t n c Q jc 



+ a T c $ c + a m 2 <j> + a Y P k $ k 



(22) 



n n c n d D icd + a - ■" <!>;. //' I), ,, - a cri$ 



where (j\ is another free parameter that controls the 
damping of the first order constraints (|15H16p . Addition- 
ally, we have introduced the familiar lapse a and shift 
(3 l functions for convenience. These are defined through 
the relations f} 1 = 7* 3 (3j (7^ being the spatial projection 



of g a b satisfying -y l lij 



-j j 9ti 



Pi) and g 1 



Equations (|17ti22[l constitute the final set of equations 
being integrated in our implementation. 



which must be satisfied for a consistent solution. An 
analogous situation is encountered with the conditions 



Cijab — ^iDjab 9jDi ab 



Cij = <)A\, 







(16) 



resulting from the commutativity of the second spa- 
tial derivatives. Henceforth we will refer to the 
{Ci, Cij, Ciab, Cijab} quantities as the C-constraints. At 
this point the resulting first order system is described 
by the evolution equations for the array of fields 
{gab, Qab, Di a b, 4>, II, together with the Z-constraints 
© and C-constraints (|15I16|) . 



D. Analysis quantities 

One way to identify the position of the boson star, 
which in the fully dynamical might prove difficult, is to 
use the energy density p = n a n h T a b- We have used the 
maximum of p, at each temporal slice to represent the 
center of each star-like configuration. 

Another useful quantity, which is of special interest in 
the case of the boson-antiboson collision, is the Noether 
charge of the theory. Due to the U(l) symmetry of the 
stress-energy tensor, there is a conserved quantity which 
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can be computed by 



N 



J°dx 6 



(23) 

We employ the Noether density J° to give additional 
evidence as to the location of the boson star. As discussed 
in [H, this quantity can be associated with the number 
of bosonic particles. 

The mass of the stars is computed by means of the 
ADM mass, which is defined as: 



M = -J- lim ( [d j9ik - d k9ij ] N k dS 

107T r— >oo J 



(24) 



where Af k stands here for the unit outward normal to the 
sphere. 

The gravitational radiation is described asymptotically 
-when the correct tetrad and coordinates are adopted-, 
by the Newman-Penrose (spin-weighted) ^4 scalar. This 
quantity is defined as 



^4 = Cabcd n' 



m n 



(25) 



with Cabcd being the Weyl tensor and where n a denotes 
the incoming null normal to the extraction worldtube 
r while m a is the complex conjugate of the null vector 
tangent to T at a constant time slice. 

To analyze the structure of the radiated waveforms it is 
convenient to decompose the signal in —2 spin weighted 
spherical harmonics as, 



^4 



l Y,, 



(26) 



where the factor r is included to better capture the 1/r 
leading order behavior of ^4. For head-on collisions of 
non-spinning objects, one could take advantage of the 
natural axis of symmetry defined by the line joining the 
centers of the objects. In such a case, the (spin weighted) 
spherical harmonics s Yi :Tn can be defined such that the 
radiation produced in the axisymmetric configuration 
would display essentially only m = modes. However, 
as we are interested in more general cases also, we adopt 
spherical harmonics defined with respect to an axis that 
can be regarded as orthogonal with respect to the orbital 
plane, which in our present case is given by the z axis. In 
this case, the radiation extracted would have non-trivial 
components for m =^ 0. 

We also focus in integral quantities that are indepen- 
dent of the specific basis of the spherical harmonics, such 
as the radiated energy. In terms of ^4, we can write 



dE 
~dt 



r 2 

ext 

Air 



* 4 (t') dt' I * 4 (i') dt' dn 

) 

(27) 

where ^4 denotes the complex conjugate of ^4, and r ex t 
is the extraction radius. Using both the decomposition 



(|26p and the orthonormalization of the spherical harmon- 
ics, this expression can also be written as 



dE 



1 



\ " 



dt 4?r ^ 

l,m 



Dt 



D, 



Q, m (t')dt'. (28) 



Finally, we notice that we will assume the resulting co- 
ordinates in the neighborhood of the extraction surfaces 
satisfy, to a reasonable level, conditions ensuring ^4 is es- 
sentially free of spurious effects (aside from those arising 
from the extraction at finite distances) [1^, [3(| . Detailed 
studies of the possible influences of these in more general 
settings will be presented elsewhere [3l|. 



III. IMPLEMENTATION ISSUES 



The first order symmetric hyperbolic system (|17ti22|) is 
implemented following techniques devised to take advan- 
tage of several useful theorems which guarantee the stable 
implementation of linear hyperbolic systems. In general 
these techniques involve: (i) discrete derivative operators 
satisfying summation by parts (SBP), which allow defin- 
ing a semidiscrete norm and an energy estimate for the 
solution [H, HE] , (ii) a method of lines employing a third 
order Runge-Kutta operator for the time integration that 
preserves the energy estimate in the fully discrete case 
[2 41 ] , and (iii) maximally dissipative boundary conditions 
that keep this energy bounded in time [13] in presence of 
boundaries. The combination of techniques (i— iii) guar- 
antee that linear problems are implemented stably and 
thus provides a direct route to a robust implementation of 
generic first order hyperbolic systems. Additionally, we 
introduce a Kreiss-Oliger type dissipative operator (suit- 
ably modified at/near boundary points so as not to spoil 
the energy estimate [26|, H3])- The introduction of dissi- 
pation allows one to control the high frequency modes of 
the solution which are always poorly described by a finite 
difference approach. In what follows we provide further 
details of the particular application of these techniques 
in the present work. 



A. Discrete operators and time evolution 

We have adopted second order derivative operators sat- 
isfying SBP and defined a discrete energy essentially as 
the L2 norm of the solution, that is 



(29) 



A I.J.K 



where XJ a = {g a b, Qab, D ia b, 4>, II, $i} is an array with 
all the evolved variables and Ua(I, J, K) is the value of 
Ua at the grid point (I,J,K). The energy estimate is 
obtained when the discrete operators reduce to the stan- 
dard three-point centered derivative operator at interior 
points while the standard (first order) two-point side- 
ways stencil at boundary points . The particular form 
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for the Kreiss-Oliger fourth order dissipation operator is 
described in [2(| [27j as is the third order Runge-Kutta 
scheme employed to advance the solution in time. 

We additionally employ the norm of the Z-constraints, 

||^||=^||^|| (30) 

to monitor the deviation of the numerical solution from 
the physical one. 

B. Characteristic structure 

As mentioned, we employ maximally dissipative 
boundary conditions. These conditions rely on knowing 
the characteristic mode structure of the principal part 
of the system. In essence, these conditions bound the 
amount of "energy" the incoming modes introduce at the 
boundaries to remain below what the outgoing modes 
carry away from the computational domain. Schemati- 
cally, this reduces to adopting conditions satisfying w~ — 
Lw + with w~ , w + denoting incoming/outgoing modes at 
a given boundary, and L a bounded matrix [23j | . As a re- 
sult, the norm (|29j) of the solution remains bounded in 
the linear regime, which is key to ensure the well posednes 
of the IB VP (for details see [II, [H, Hi] ) . In our particular 
case, we choose L = thus setting the incoming modes 
to zero. The characteristic decomposition (eigenmodes w 
with velocity v) for the EKG system (|17H22|) is given by 

w ab ° = g ab , v° = (31) 

Wl J = D iab - N k D kab Mi , = -A4 I3 k 
w ab (±) ee Q ab - era g ab ± Af k D kab , v ± = -A4 P k ± a 

with Af k the outgoing normal vector to the boundary sur- 
face and the supra-indices {0,/3, ±} denote modes prop- 
agating with the different speeds {v°,v l3 ,v ± }. Thus the 
incoming modes are given by u>~, and Wi a \P as long as 
A4 P k < 0. To apply the maximally dissipative condi- 
tions at the numerical level, we follow the prescription 
described in [12, HI] and simply enforce 

dtwab^ = 0, (32) 
dt^L = 0, if A4 (3 k < 0. (33) 

C. Mesh Refinement 

Adaptive mesh refinement (AMR) provides increased 
resolution where and when needed, and therefore be- 
comes an essential tool in resolving a wide range of dy- 
namical scales, both spatial and temporal. As described 
by Berger and Oliger [321 ] , a given numerical grid can be 
further refined based on dynamical measures of the nu- 
merical error by the creation of refined grids which are 
evolved in sync with the coarse grid. We have built our 



code taking advantage of the computational infrastruc- 
ture called HAD (an outgrowth of the code presented in 
[liJHij])- HAD implements a (slight) modification of the 
standard Berger-Oliger strategy that guarantees preserv- 
ing the stability properties of the unigrid implementation 
and significantly reduces spurious reflections at interface 
boundaries [35] . In HAD creation or destruction of finer 
grids is automated so that points which display an er- 
ror above some threshold are "covered" by a finer grid. 
The error associated with each point is computed via a 
self-shadow hierarchy [36| . In this technique, the error is 
computed on any given grid by subtracting the solution 
on that grid with that on the next coarser level. Because 
these solutions evolve with different resolutions, their dif- 
ference represents a measure of the local truncation error 
without recourse to more complicated shadow schemes. 
In particular, at any given point, the L2 norm of the 
difference between the solutions over all evolved fields 
is taken as the truncation error estimate. Although we 
have performed binary boson star runs employing AMR, 
in this work we take only partial advantage of the mesh 
refinement capabilities of HAD. One reason for this is 
that we want to compare significantly different scenarios 
and we find it useful to ensure similar discretizations are 
employed across the cases considered. Our procedure is 
to let the infrastructure define the grid structures at the 
initial time and keep it fixed throughout the simulations. 
The truncation error value employed by the self-shadow 
hierarchy to define the grids is chosen such that the origi- 
nal stars are contained within the finest grid and that this 
grid covers the whole region in between. Fixing the grid 
structure at the initial time reduces some of the overhead 
and memory required by dynamic regridding. 



IV. SIMULATIONS AND RESULTS 

In this section we describe the simulations we have 
performed both to check the correctness of the imple- 
mentation and investigate interesting aspects of boson 
star head-on collisions. We begin in llV Ai by considering 
a single isolated boson star configuration and verify the 
code is able to reproduce known features of the solution. 
The evolution of the star is performed employing two dif- 
ferent coordinate systems, a static one and the harmonic 
one. The former coordinate system is the natural one 
employed in constructing a single static boson star solu- 
tion, and so the evolution with it should be trivial. The 
latter one does not conform to the ansatz employed to 
obtain the solution and some non-trivial dynamics can 
be expected. 

In Section IIVBI we consider different binary systems 
described by boson stars with zero initial momentum. 
Initial data for these systems is obtained by a simple 
superposition of a single boson star with another which 
is identical to it up to a possible phase or reflection of its 
natural frequency. The stars are placed sufficiently far 
from each other to ensure the constraints are satisfied to 



6 



a degree consistent with the value obtained in the single 
boson star case. 

We have performed three different types of head-on col- 
lisions: collisions of equal boson stars IIVB 11 collisions of 
boson stars whose phase are in opposition IIV B 21 and fi- 
nally collisions of boson stars with antiboson stars If V B 31 
(these differences are explained in detail below). 

For all of the simulations presented here we have 
adopted damping parameters do = a \ — 1- a Courant 
factor At /Ax — 0.25 and Kreiss-Oliger dissipation pa- 
rameterized by 0.03 (as described in [2^, [23]). The sim- 
ulations have been performed on 64 to f28 processors 
employing 3 refinement levels, placed as defined by the 
shadow hierarchy. 



A. The single boson star 

As described in the introduction, a boson star is a self 
gravitating configuration of scalar field with Lagrangian 
given by Jl}. In a particular coordinate system, one as- 
sumes a particular form for the spherically symmetric 
scalar field 



= 4>o(r) e 



(34) 



Therefore the field values at any given point oscillate with 
a constant frequency u>. Details of the calculation of 4>o( r ) 
are presented in Appendix[XJ The geometry, on the other 
hand, remains static in these same coordinates. The re- 
sulting spacetime provides a reasonable model for that 
resulting from a compact object in equilibrium. 

We have carried out simulations of single boson stars 
under two coordinate choices. In IIV A II we summarize 
the results in the coordinate system on which the solu- 
tion is static and the results in the harmonic gauge is 
presented in IIV A 21 The simulations presented for both 
cases employ a grid structure consisting of three levels of 
refinement. The finest level covers the radius R95 of the 
star (defined as region of space containing approximately 
95% of the star's mass) and the other boxes surrounding 
it with half the resolution each respectively. For all these 
runs the star is located at x % — (i — 1..3) and the 
computational domain is given by x l € [—144, 144]. 

To obtain a measure of the convergence rate of the so- 
lution, we have performed simulations with different res- 
olutions for the boxes. The finest grid of each resolution 
considered has Ax £ {0.375,0.25,0.156} for the static 
coordinates and Ax € {0.5, 0.375, 0.25} for the harmonic 
ones. The others grids, as mentioned, have a gradual 2:1 
ratio in their discretization length. 

Below we detail certain convergence tests of the code 
about which we would like to make a general comment. 
Residuals of constraints are expected to approach zero 
with increases in resolution while field variables should 
approach a unique solution, and this behavior is indeed 
what we find. The rate of convergence can also be ex- 
amined and compared to the expected approximation or- 
der of the difference equations solved, and we found the 
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FIG. 1: Single boson. Phase oscillation of the real part of 
the scalar field at the center of symmetry up to t — 50 for 
the resolution Aa; = 0.25 (similar behavior was followed up 
to t — 500). The solid line indicates the analytically expected 
value, the crosses show the values read-off from the numerical 
solution within the static coordinates and the circles show the 
same by using the harmonic coordinates. 



expected second order convergence in the unigrid case. 
However when dealing with a complicated grid hierar- 
chy, precise measurements of this rate become somewhat 
problematic, especially when the hierarchies produced 
vary with the truncation error threshold chosen. Also, 
we should mention that more detailed convergence tests 
with this code infrastructure applied to other problems 
demonstrate the expected order of convergence [H, H3 ■ 



1. Single boson star in static coordinates 

In order to maintain the coordinates such that the 
spacetime remains static the H a (t,x l ) functions in (JT]) 
need to been chosen in a suitable way. Since we are con- 
sidering a static situation, the condition on H a can be 
read directly from the equation ([9]) by imposing the con- 
straints (i.e., Z a — 0) and the staticity of the spacetime 
(i.e., dtg a b = 0). In a explicitly static scenario, this initial 
value of the H a (t = 0, x l ) is preserved for all times. 

We begin with an initial data describing a star where 
the value of the scalar field defining it at the origin is 
\<t>a{r = 0)| = 0.01. This star, which is obtained by 
solving the initial data problem with a high resolution 
one-dimensional code as described in the Appendix [3] 
has an ADM mass of Mid = 0.361 and radius R95 = 
19.6, and lies well inside the stable branch of solutions. 

The evolution of an isolated star displays the expected 
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FIG. 2: Single boson. Loo of g xx versus time for solution ob- 
tained employing the static coordinate condition. The figure 
displays the result of the evolution for three different base- 
resolutions Ax = {0.375,0.25,0.156} (in addition, each nu- 
merical solution is obtained with three levels of refinement). 
The values obtained converge to a constant value as expected, 
and the maximum variation observed to t = 50 is ~ 0.0073% 
for the coarsest base-resolution while for the highest one is 
just ~ 0.0015%. 



behavior. In particular, the geometry is static and the 
behavior of the scalar held agrees with Eq. (j34|) as demon- 
strated in Fig. [T] where the central value of the real part 
of the scalar field versus time describes an oscillatory 
behavior with a frequency u> = 0.96 ± 0.03. This is in ex- 
cellent agreement with the frequency obtained with the 
one-dimensional initial-data problem code uj\d = 0.976. 



2. Single boson star in harmonic coordinates 

In this section we revisit the initial data of the previ- 
ous one, but in this case we do not adopt coordinate 
conditions that explicitly demonstrate the staticity of 
the spacetime. One goal of this exercise is to test the 
ability of the harmonic-coordinate condition to adapt to 
the physical problem under consideration. This will be 
important for the binary cases considered later, as that 
problem is certainly not static. 

The implementation evolves the system for as long as 
the code was run without any signs of instabilities. For 
this particular case we have evolved the system for more 
than t = 1500. During this evolution, we have monitored 
the ADM mass which we extract at r ext = 100 and obtain 
M = 0.364 which decreases by less than 3% / 1% for the 
coarsest/finest resolution considered (Ax = 0.5 / 0.25 for 



the innermost grid) by the time we stop the code. This il- 
lustrates the ability of the code to preserve this conserved 
quantity for long times. 

We have also checked that for this particular solution, 
both the harmonic and the static coordinates agree at 
the center of symmetry. This implies that the scalar field 
should have the same local oscillatory behavior (f3~4"|) at 
r = in both coordinate systems. This is indeed the case 
and is illustrated in Fig. [TJ In particular, the measured 
oscillation frequency is ui = 0.97 ± 0.06, Additionally, we 
have checked that the absolute value of scalar field at 
the center of the star varies by at most 12% / 2% for the 
coarsest /highest resolution considered. 

Coordinate effects do arise however. As men- 
tioned above the geometry should have some non-trivial 
coordinate-induced dynamics since we are not adopting 
a coordinate system on which the solution is explicitly 
static. This effect can be seen in Fig. [3] where the evolu- 
tion of the maximum of function of coordinate 
time t for 3 different resolutions is displayed. As evident 
in the figure, there is an initial transient variation of the 
metric value which later approaches a constant value. 

Finally, we monitor the constraints in Fig. [¥] The fig- 
ure shows the L 2 norm (f3"0"|) of the physical constraints 
(i.e., the Z-constraints ©) and its behavior when the 
resolution is improved, converging clearly to zero. The 
measured convergence rate is 2.9 between the low-and- 
medium and 2.8 between the medium-and-high resolu- 
tions. 



B. Binary boson-stars head-on collisions 

Initial data for the collision of two boson stars is 
obtained by a simple superposition of the single boson 
solution in a way described below. Additionally, in order 
to consider different configurations, we take advantage 
of the following observations: 

(i) The solution of the initial data problem is unaf- 
fected by a phase difference in the ansatz assumed for 
the scalar field. Namely, (f> — (j>Q(r)e l ^ wt+9 ^ gives rise to 
the same initial data for {g a b, dtg a b = 0}. 

(ii) if {g ab , d t g a b = 0, d t cj) = iuxf> } provides 
consistent initial data, so does {g a b, 4>o, dtg a b = 0, <9t</> = 
—itucpo}. The difference is solely in a frequency-reflection 
of the boson star and is known as an anti-boson star. 

We will exploit these observations to consider what we 
will refer to as a boson in phase opposition by taking 
9 = 7T in (i) and/or an anti-boson in (ii). 

The initial data we consider is schematically repre- 
sented by considering the following construction 

6 = (1) (x - xi,y, z) + </)( 2 \x - x 2 ,y, z) (35) 
ip = i/;( 1 \x-x 1 ,y,z) + ip( 2 \x-X2,y,z)-l 
a = a*- 1 -* (x — x\, y, z) + (x — x 2 , y, z) — 1 
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FIG. 3: Single boson. L x of g xx versus time for solution 
obtained employing the harmonic coordinate condition. The 
figure displays the result of the evolution for three different 
base-resolutions Ax = {0.25,0.375,0.50} (in addition, each 
numerical solution is obtained with three levels of refinement). 
After a small transient behavior, the values obtained converge 
to a constant value. The maximum variation observed to 
t — 500 is ~ 0.72% for the coarsest base- resolution while 
~ 0.11% for the highest one. 



where denotes the corresponding field of the boson 
i, centered at (^,0,0) and the value of (jy^> will be dic- 
tated by the type of boson star considered, that is: a 
boson star, a boson star in phase opposition or an anti- 
boson star. Notice that three fields (the scalar field (f>, 
the conformal factor ip and the lapse a) are enough to 
set an initial data consistent with the EKG equations, 
as it is described in the appendix. Under this approach 
we are making the assumption that the boson stars are 
described with a single global scalar field cf) instead of 
considering two sets of complex scalar fields to represent 
each star. This choice fixes in a straightforward way the 
interaction of the stars. 

In order to consistently choose the initial position of 
the boson stars we have measured, at the initial time, the 
Lqo norm of the Hamiltonian constraint for different coor- 
dinate separations of the centers D. Since at sufficiently 
large distances the adopted data satisfies the constraints 
by construction, the constraint's behavior versus D pro- 
vides a measure of the distance that must be adopted. 
As illustrated in Fig. the error decreases very rapidly 
with separation and we have chosen the value D = 50 
(which corresponds to x\ — —X2 — 25) that is nearing 
the asymptotic behavior. In addition, we will use for all 
the head-on collisions the same form of the scalar func- 
tion 4>o{r) in order to compare easily the different cases. 



FIG. 4: Single boson. Convergence of the L2 norm of the Z- 
constraints residual (defined by equation (|30p ) versus time by 
using harmonic coordinates, for the same resolutions than in 
Fig. [3] As the resolution is increased the constraint violation 
is clearly reduced. 

1. Merging of the boson/boson pair 

In this case, two identical boson star configurations 
are adopted to define the initial data. Since the stars 
adopted have no initial velocity, their initial behavior is 
marked by a slow approach towards each other as they 
feel their gravitational attraction. As the evolution pro- 
gresses however, the resulting behavior depends on the 
mass of the single initial boson stars. Previous studies 
which concentrated on a particular set of initial masses, 
always displayed a collision of stars giving rise to the 
formation of a black hole [l2T |. More recently, the work 
of [1, [38| , which considered boosted stars with large ki- 
netic energy, presented cases where the stars appear to 
pass through each other in a sort of "solitonic" behavior. 

Here we want to study the situation in which there is 
a final regular object (that is, containing no singulari- 
ties). To do so we begin with a broad parameter search 
along the masses of the boson stars searching for cases 
where the final compact object does not collapse to a 
black hole. As shown in Fig. [6l for stars with masses 
around M — 0.26 (corresponding to an amplitude of 
4>o(r = 0) = 0.005) the final object appears to avoid 
black hole formation. This avoidance is indicated by the 
geometric variables tending to a smooth slow oscillatory 
behavior for M < 0.26 while a marked violent behavior 
is displayed for larger values, together with a collapse of 
the lapse function. 

We thus focus on the particular case where each star 
has masses M = 0.26, which is approximately 40% of 
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FIG. 5: Boson/Boson pair. Loo norm of the Hamiltonian 
constraint residual versus the initial (coordinate) separation 
of the boson stars. As the distance is increased the violation 
decreases as expected since each star is defined so as to satisfy 
the constraints. We choose a value of D — 50 which lies close 
to the asymptotic behavior and is still manageable in terms 
of the boundaries. 



the maximum allowed mass on the stable branch for the 
potential (fl2|) . The radius of the single star is around 
i?95 = 27, and in this case we extend the computational 
domain to x % E [—320, 320]. As mentioned above, the re- 
finement regions are rectangular boxes covering the cen- 
ters of the stars and the span between them with a grid 
spacing at the region of the star of Ax = 0.50 for the 
lowest resolution and Ax — 0.375 for the highest one. 
We compute the values of the ADM mass and "J" 4 at ex- 
traction surfaces located at r ext — 140, 170 and 200, 
where the grid spacing is given by Ax = 4 for the lowest 
resolution and Ax = 3 for the highest one (i.e., there 
are 3 levels of refinement). Since the wavelength of the 
observed radiation is of the order A ~ 100, the travel- 
ing wave is well represented by the grid-structure at the 
extraction location. 

In an effort to ensure fidelity to the continuum prob- 
lem, we monitor constraint residuals and convergence of 
the metric variables for these evolutions. Notice however 
that, again, the grid structure for different resolutions 
differs and so our evolutions do not lend themselves to a 
traditional convergence study. However we do see the ex- 
pected behavior, in that the constraint residuals decrease 
in regions where the resolutions differ and the fields indi- 
cate convergence to a common value. This is illustrated 
in Fig. [7] where the maximum value of g xx versus time is 
displayed. As evident in the plots, the values for different 
resolutions agree quite well until a time of about t ~ 800 



FIG. 6: Boson/Boson pair. Geometry behavior for differ- 
ent values of central value of the scalar field that defines the 
stars. As the stars merge, the behavior of geometric vari- 
ables changes drastically for different values. Stars with ini- 
tial central densities > 0.005 seem to give rise to a black hole 
as illustrated by their considerable growth; those below this 
value yield much smoother behavior. This is illustrated in the 
figure which displays the maximum of the component of the 
metric g xx versus time for different central amplitudes of the 
scalar field {<j> = 0.02,0.015,0.01,0.05} respectively. These 
give rise to stars with to masses {M = 0.47, 0.42, 0.36, 0.26} 
respectively. 



where accumulation of errors and boundary effects be- 
come apparent. 

To illustrate the dynamics displayed by the solution, 
we present in Fig. [H] and a sequence of plots illustrating 
the behavior of the Noether density J°. The extremes of 
this function correspond, at least initially, to the centers 
of the stars and as the evolution proceeds, the extremes 
give an indication as to the movement of the stars. Fig. [5] 
presents two-dimensional slices (at z = 0) of J°. Notice 
that the maximum value of J° is reached after the col- 
lision. This suggests that the boson stars have merged 
and oscillations of the merged object are clearly distin- 
guishable. Fig.[5]shows the nature of these oscillations in 
more detail by presenting contour plots of</° = 2xlO~ D 
at more frequent times. 

The gravitational radiation produced in the collision, 
as encoded in $4, is presented in FiglTUl which illustrates 
the (real) coefficients of its Z = 2 spin-weighted spherical 
harmonic modes. Due to the symmetry of this partic- 
ular problem and with the spherical harmonics defined 
with respect to the z axis, there are the following simple 
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FIG. 7: Boson/Boson pair. Loo(g xx ) as function of time for 
three different resolutions and their differences. We can ob- 
serve clearly that the merger occurs at approximately t ~ 550. 
The solution is qualitatively convergent during the first part 
of the merger until t ~ 800 when boundary effects negatively 
affect the convergence. 



relations between the non-trivial coefficients C 2 



Re{C 2>2 } 
Re{C 2:2 } 



Re{C 2 ^ 2 } 
-^/3/2Re{C 2fi } 



(36) 



Numerically we have found that Re{C 2 . 2 ] = Re{C 2y - 2 } 
and Re{C 2y2 } = —1.22 Re{C 2y o}, which is in excel- 
lent agreement with the expected analytical relations 
The gravitational radiation emitted during the collision 
reaches the first extraction surface (at r ext — 140) at 
t = 600. It takes around t = 200 — 300 to reach the 
boundary, be reflected and pass again through the same 
extraction surface, as can be seen in Fig.QT] Notice that, 
despite having the surface extraction far from the sources, 
they are still not well within the "wave-zone" since r\T/4 
displays still a slight dependence on r ext . Nevertheless, 
the structure of the radiated wave is evident in the plot. 

In Fig. [T2] we have plotted the emitted energy (or 
power) dE/dt (eqn. [28)) as a function of time for dif- 
ferent resolutions. We find that practically all energy is 
radiated in the I = 2 mode which would simplify the pa- 
rameterization of the obtained waveforms for their use in 
data analysis searches [30(. 

These evolutions therefore suggest that the merger of 
the two boson stars produces an oscillating single boson 
star. However, we note that the resulting star does not 
correspond to the boson star with mass equal to the sum 
of the individual masses. One might suspect the collision 
dispersed energy in the form of scalar and gravitational 




t = 300 




t = 900 




FIG. 8: Boson/Boson pair. 2D 2 = cuts of the Noether 
density J° at different times. As the stars come closer and 
merge, the maximum value of J° grows significantly followed 
by a quadrupolar oscillation. 



radiation, but that does not appear to be the case. The 
ADM mass hardly decreases and the gravitational wave 
output is minimal. Instead, it would appear that the 
resultant object has yet to settle into a stationary star. 



2. Collision of boson/boson in opposition of phase pair 

The previous case considered, in its early stages, can 
be regarded as generic for equal mass objects in a 
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FIG. 9: Boson/Boson pair. Contours corresponding to the 
value J° = 2 x 10 -5 at the z = plane. From top to 
bottom the contours are shown for times {0,160,240,300}, 
{340, 380, 420, 460},{500, 540, 580, 640}, {680, 760, 840, 900} 
and {980, 1060, 1140, 1200} indicated by solid, solid-with- 
crosses, solid-with-squares and solid-with-triangles respec- 
tively in each plot. As the stars get closer, the initially spher- 
ical contours deform until merging as illustrated by a cusp 
on the top-most figure. Afterwards, J° exhibits essentially 
quadrupolar-type oscillations. 

head-on configuration. The inner-structure details only 
become relevant when the objects become sufficiently 
close to each other as indicated by the effacement the- 
orem [nm. 

However, as the stars approach each other, the par- 
ticular details of the star under consideration can have 
strong consequences. Let us consider the more general 
possibility of having a relative phase difference between 
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FIG. 10: Boson/Boson pair. Coefficients corresponding to the 
I — 2 modes of r^4 as a function of time, extracted at r ex t = 
200. After some initial transient due to spurious radiation in 
the initial data, the signal is clearly visible corresponding to 
the merger of the stars and later decaying to smaller values. 
At late times, t > 800, contamination with boundary effects 
obscures the extracted signals. 

the single boson stars that we will use to construct the 
global solution (|35j). as for instance, 

0« = Mr) (37) 
0( 2 ) = 0o ( r ) e -*(«*+«) (38) 

with 9 the relative phase. We will concentrate here on 
the extreme case, 9 = tt and so the stars are in phase- 
opposition. Recall that while Einstein equations are in- 
sensitive to this phase difference, the Klein-Gordon equa- 
tions are. Notice that in our head-on configuration the 
surface x — is an anti-symmetry plane for the scalar 
field, so it must remain zero in that plane. In the present 
case, the observed behavior can be regarded, to a certain 
extent, as being influenced by a "repulsive" interaction 
of two objects trapped in a gravitational potential well. 
This leads to the objects oscillating around the coordi- 
nate position {x = ±15, y = 0, z = 0). Elucidating the 
final fate of this system will require much longer evo- 
lutions (making sure there is no causal interaction with 
the boundaries). Additional comments about this case in 
light of result obtained in the other ones will be discussed 
in the conclusion section. 

The gravitational radiation produced in this scenario 
is considerably weaker than the previous case, as the in- 
volved objects never acquire significant velocities to in- 
duce a strong-time dependent quadrupole. In fact, our 
extracted radiation is of the same order as the spurious 
signal produced by the initial data. From these results 
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FIG. 11: Boson/Boson pair. The coefficient C2,+2 for dif- 
ferent extraction radii. This figure illustrates both the out- 
going waves due to the dynamics of the spacetime together 
with incoming waves that have traveled to the boundary and 
bounced off it. Additionally some remnant dependence on r 
is visible indicating the extraction is still performed not suffi- 
ciently far from the sources. Nevertheless the structure of the 
outgoing waves is clearly visible. 



FIG. 12: Boson/Boson pair. Radiated energy as a func- 
tion of time. The top plot indicates the energy radiated in 
all modes (dashed line) together with that radiated solely in 
I — 2 modes. Clearly, until boundary effects begin to affect 
the results, practically all energy is radiated in the I — 2 
modes. The bottom plot displays the difference between the 
I = 2 radiated energies at three different resolutions, indicat- 
ing convergent behavior until t ~ 950. 



coupled with the analysis presented in [El we conjecture 
that the maximum of the radiation happens at 9 = 
while the minimum corresponds to 9 = tt being the radi- 
ation a smooth function of the phase difference. However 
this must be substantiated by studying several represen- 
tative cases. 



3. Boson against antiboson 

Another interesting case is the collision of a boson star 
with an otherwise identical star except with the oppo- 
site charge density. Such a star is called an antiboson 
star and rotates in the opposite direction as its partner 
in the complex plane. Recall that the initial value prob- 
lem solution was degenerate upon the reflection lu — > —lj. 
Additionally, while this change leaves the geometry un- 
changed, the associated Noether change is affected by an 
overall sign change. In this section, we study the dynam- 
ical behavior of a binary system composed of a boson 
and an antiboson star. The initial conditions for such a 
scenario are simply obtained with the following choice, 

4F> = Mr) e~ 1 ^ (39) 

^ = Mr) e +luJt ■ (40) 

We obtain the evolution of such a system, and as with 
the phase-opposition case, the early behavior agrees with 



that of the boson-boson case. As time progresses how- 
ever, notable differences arise which are illustrated in 
Fig. [THl The dynamical behavior when the stars are close 
is less strong than in the boson-boson case but more 
so than in the boson-phase opposition case. As a re- 
sult, there is non-trivial generation of gravitational waves 
though their time-scale is longer than and their strength 
smaller than the boson-boson case. The total radiated 
energy however, is similar as illustrated in Fig. 1151 



V. CONCLUSIONS 

We have studied the behavior of boson stars both in 
isolation and head-on collisions with a three-dimensional 
implementation of the Einstein equations. Studies of the 
isolated star in Section |IV3 provided a non-trivial test of 
the implementation and showed it capable of accurately 
extracting delicate features of the solution such as the 
oscillation frequency. Additionally by considering a co- 
ordinate condition not adapted to the static spacetime 
under study, the solution obtained further illustrated the 
ability of the harmonic coordinates to deal with compact 
objects in a smooth manner. 

We then considered head-on collisions of different con- 
figurations of boson stars. Our studies revealed the inter- 
esting and varied phenomenology interacting boson stars 
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FIG. 13: Boson/Boson in opposition of phase pair. 2D z — 
cuts of the Noether density J° at different times Although 
the stars come closer, they never seem to merge. 
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FIG. 14: Boson/Boson in opposition of phase pair. Con- 
tours with value J° = 2 x 10 -5 on the 2 = place. 
From top to bottom the contours are shown for times 
{180,360,540,720} and {720,900,1080,1200} indicated by 
solid, solid-with-crosses, solid-with-squares and solid-with- 
triangles respectively in each plot. The stars initially get 
closer (top frame) but afterwards they move apart from each 
other (bottom frame), this process continues on as the stars 
are trapped in a common gravitational well. 
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FIG. 15: Boson/ antiBoson pair. The C2,m coefficients of the 

can give rise to. Additionally, these scenarios provided r* 4 decomposition (top frame) and the emitted radiation ver- 

an excellent test of the implementation under strongly sus time (bottom frame) . 
dynamical conditions. This implementation is being em- 
ployed to deal with black hole spacetimes, and the results 
will be communicated elsewhere. 

In addition to the notable differences obtained, it is in- proper distance between the origin and these maximums 

teresting to pay a closer look at the position of the boson (i.e., that can be identified with the center of the boson 

stars as a function of time. The positions are determined stars) located on the +x direction. Due to the symmetry 

by the location where the energy density reaches a max- of the problem this is enough to draw conclusion on both 

imum in a given neighborhood. Fig. [T5] shows the max- stars. Included in the figure is the position a star would 

imums of the energy density p while Fig. [TO] shows the have as dictated by a simple Newtonian behavior of two 
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FIG. 17: Boson/antiBoson pair. Contours with value 
J° = 2 x 10 -5 on the z = place. From top to bot- 
tom contours are shown for times {180,360,540,600,660} 
and {660, 780, 840, 900} indicated by solid, solid- with-crosses, 
solid-with-squares and solid-with-triangles respectively in 
each plot. Only the right star has such a contour initially 
as the left has a negative initial value for J°. As the stars 
come closer to each other, the "tunneling" behavior is illus- 
trated by the positive value appearing on the left side (top 
frame). As time progresses, this process essentially reiterates 
itself from the left side to the right side (bottom frame). The 
dynamics that follows repeats this side to side 'motion' in the 
Noether density with the stars trapped in a common gravita- 
tional well. 



t = 900 




t= 1200 




FIG. 16: Boson/antiBoson pair. 2D z = cuts of the 
Noether density J° at different times Although the stars come 
closer and merge, the bosonic part is distinguishable from the 
antibosonic one by means of the different sign of the Noether 
density. 



point-particles with the mass of the boson star. 

A priori, since the initial separation of the stars is 
of the same order of their radii, non-Newtonian effects 
would be expected to become relevant early on. This 
is clearly visible for the boson-boson and boson-boson 
in phase opposition cases, however, the Newtonian re- 
sult is a reasonably good approximation for the boson- 
antiboson collision. This results from the significantly 
different interaction between their scalar fields in these 
cases, while in the boson-boson case a merger takes place 
at relatively early times, in the boson-boson in phase op- 
position the stars approach and essentially go through 
each-other. 

A further remark about the particular cases that we 
have evolved is that the difference between the coordi- 
nate and the proper distances, during most part of the 
evolution (the early part), is given in all cases by a slowly 
varying function of time, compared to the time scale of 
the solution. 

Additionally, these plots also explain, qualitatively, 
why the radiated energy is stronger for the boson-boson 
case since it is the one where the source's quadrupole 
has the most dynamic behavior. It would be tempting 
to regard this as further evidence that in rather simple 
scenarios, the coordinate behavior provides a reasonable 
description of the objects dynamics as seen in the con- 
text of equal-mass binary black hole cases jU [42], HH, |44| . 
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FIG. 18: The maximum of the energy density p is plotted as 
a function of time for the three different cases. 



However one must be cautious in the present case as we 
are considering a particular model to define the scalar 
fields through a single complex field. 

Also, interesting differences are observed when com- 
paring our results with those of the study of Q-balls [3] • 
Here, the inclusion of gravitational effects modify the dy- 
namics strongly. Besides not needing to boost the stars 
for a collision to take place, the gravitational attractions 
of the stars in our case is evidenced by the "trapping" 
of the stars instead of "bouncing apart" (for the boson- 
boson in phase opposition) or going through each other 
(for the "boson-antiboson" ) case as in [l4[ ■ 

Finally, the results presented here illustrate the rich 
behavior displayed by binary boson star dynamics even 
in the rather simple scenario of a head-on collision. In 
particular, the dependence on possible phase differences 
revealed these phases strongly influence the outcome of 
the problem. In appendix [B] we present a simple treat- 
ment that sheds light in the influence different phases or 
frequency-reflection in one of the stars' description have 
in the dynamics of the problem. In particular, this treat- 
ment suggests that energy arguments explain why the 
boson-boson case gives rise to a merge while in the boson- 
opposition in phase boson case they seem to repel each 
other. 

Having examined the landscape of head-on boson star 
collisions behavior, near future studies will concentrate 
on the more complex situations of orbiting cases as well 
as boson-star black hole binaries. 



FIG. 19: The proper distance from the center of the boson 
star to the center of the computational domain as a function of 
time for the different cases and the Newtonian approximation. 
The position of the boson is identified with the maximum of 
the energy density. 
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APPENDIX A: INITIAL DATA FOR AN 
ISOLATED BOSON STAR 

The initial data for the boson star configuration is com- 
puted first in spherical symmetry with a one-dimensional 
code. The resulting solution is then extended to three 
dimensions. The one-dimensional solution is obtained in 
the following way. First, we adopt a specific ansatz for 
the scalar field, 

^(t,r)=(t> (r)e(- iut l (Al) 

With this assumption, our goal is then to to find <p Q (r) , U) 
and the metric coefficients such that the spacetime gen- 
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erated by this matter configuration is static. We begin 
by considering the problem in polar-arcal coordinates as 
is done in, for instance, @, H, 0]. The line element in 
these coordinates takes the form, 



ds 2 = -a (rf dt 2 + a {r) 2 + r 2 dn 2 . 



(A2) 



The equilibrium equations in this coordinate system are 
then given by, 



a 2 - I 
r 

-\-4irr 

a 2 - I 

r 

4- 47rr 



,2i2 



+ m a 



m 2 a 2 6 2 



< = - (1 + a 2 - 47rrVm 2 ^) ^£ 
u 2 



a' 



m ) 4> a 



(A3) 

(A4) 
(A5) 

(A6) 



where a prime denotes differentiation with respect to r. 
In order to obtain a solution of this system adapted to 
the desired physical situation we provide the following 
boundary conditions: 



o(0) 


= 1, 


(A7) 


0o (0) 


= 0c, 


(A8) 


t>o (0) 


= o, 


(A9) 


d r a (0) 


= o, 


(A10) 


lim a (r) 

r — >oo 


= 1, 


(All) 


lim (j) (r) 

r — >oo 


= o, 


(A12) 


lim a (r) 


= i; 


(A13) 



which guarantee regularity at the origin and asymptotic 
flatness. For a given value of central value of the field 
(j) c these equations and boundary conditions only ad- 
mit solutions for a discrete set of values of w. In our 
particular case, we we are interested on the fundamen- 
tal (lowest frequency) solution. The problem is solved 
by integrating from r = outwards using a second or- 
der shooting method implemented using the numerical 
package LS0DaJ45|. The equations are integrated setting 
a(r = 0) = 1 initially, after all the equations are solved, 
the lapse is rescaled in order to obtain a function which 
asymptotes to 1 at infinity [48|. The same rescaling is 
then performed to the frequency to. 

Once the solution is computed in this coordinate sys- 
tem a change of coordinates is performed to maximal 
isotropic ones: 



In these coordinates the extension to three dimensions is 
direct since the solution is explicitly spatially conformally 
flat. We perform the transformation of coordinates as in 
Q. Since air) and ip (f) (and the scalar field <p (f)) 
are scalar functions of the spacelike hypersurface, and 
the time-slicing is not changed, the extension to three- 
dimensions involves only writing these functions as func- 
tions of x, y and z such that f 2 = x 2 + y 2 + z 2 . This task 
is performed using a 5 point Lagrange interpolation. 

This way we obtain initial data for g a b and <f>, the rest 
of the fields for the three dimensional code are chosen 
as follows: Q a {, = 0, II is computed from ansatz (|A1[) 
and Di a b and <&i are calculated using constraint equa- 
tions (|T4"|) . That completely defines the initial data for a 
boson star. 



APPENDIX B: ENERGY CONSIDERATIONS 

To gain some physical insight into these results, we 
consider here how the energy density behaves for the 
three different cases studied. To simplify the discussion, 
we consider the energy density in flat spacetime, i.e. we 
assume a Minkowski metric. Under this assumption, the 
energy density for a complex scalar field 4> can be written 
as 



P 



mi 



(Bl) 



where II and $ are the derivatives of <fi defined as before 
in Eqs. (fT3"|) and (|14p . Now, we treat the scalar field as 
given by the superposition 4> — 4>i + (j)2 with 4>i,4>2 de- 
scribing the different stars we study. The energy density 
associated with this superposition can be expressed as 



p = pi + P2 + A 



(B2) 



ds 2 



-2 (~ 



(f) dt 2 + V> 4 (f) (dr 2 + f 2 dft 2 ) 



where pi the energy density that corresponds to the field 
<f>i in isolation and A is the interaction potential which 
vanishes when the stars are well separated. The interac- 
tion potential can then be written as 

A = - [nii±2 + nif± 2 + $i$2 + $1*2 

+m 2 (<M 2 + </>i0 2 )] . (B3) 

We next assume the scalar field assumes the form 

6, = 0S(t,a; < )e iu * (B4) 
2 = 0°(i,:rV (eart+9) , e = ±l (B5) 

where <j)\ represents the complex field configuration asso- 
ciated with a dynamic boson star and 2 represents the 
field configuration of the other star. Both stars are char- 
acterized by the same frequency ui but the second star can 
be a regular boson (e = 1, 6 = 0), a boson in opposition 
of phase (e = 1, 6 = tt), or an antiboson (e = — 1, 9 = 0). 
The profile of each is described by the real field 0^(t, x l ) 
(A14) which is a assumed to be a slowly varying function of 
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time so the most dynamical part of the time dependence 
is represented by the harmonic factor e luJt . 

Under these assumptions, a straightforward computa- 
tion reveals that the effective interaction potential can 
be expressed as 

A = A cos [(1 -e)wt-6] . (B6) 

Notice that when the two stars are centered in the same 
position the function Ao is strictly non-negative, since it 
is a sum of quadratic terms, 

A - [n°n° + n°n° (B7) 

+$?$° + n<t° 2 +m 2 + , 

and the behavior of the interaction term for each case 
considered is governed by that of cos [(1 — e) u)t — 9]. 
Thus, for the problems analyzed in the present work, the 
interaction potential for each case, is, 

A B -b = A , (B8) 
Ab-opb = — A , 
Ab-ab = A cos(2wi) . 



Since the contribution of Ao is positive for the boson- 
boson case, the resulting gravitational potential associ- 
ated with this case would be deeper than if not present. 
As a result, a merge of the two stars would seem a natural 
consequence. On the other hand the Ao contribution has 
the opposite sign for the boson-opposition phase boson 
case, translating into a gravitational potential exhibit- 
ing a bump at the center with two minimums around it. 
Such a bump would suggest it energetically preferable for 
the stars not to merge. Finally, the boson-antiboson in- 
teraction potential is governed by a varying function of 
time. If the interaction time-scale is much longer than 
~ 1/cj, the interaction term essentially integrates away 
to 0. Whether this is the case depends on the collision 
under study and, in particular, the stars' momentum as 
they come close to each other. 
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